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Abstract 


We report on progress in the development and application of a coupled Boundary 
Element/Finite Volume Method temperature-forward/flux-back algorithm developed to solve 
conjugate heat transfer arising in 3-D film-cooled turbine blades. Here, heat conduction within 
the blade is coupled to heat transfer in the external fluid flow field that is convecting heat into/or 
out of the blade. In the BEM formulation no interior mesh is generated and the surface heat flux 
is computed in the solution. We adopt a loosely coupled strategy where each set of field 
equations, Navier-Stokes for the external field and heat conduction for the internal field, is 
solved to provide boundary conditions for the other. The equations are solved in turn until 
iterative convergence criteria requiring continuity of temperature and heat flux are met at the 
fluid-solid interface. The NASA-Glenn turbomachinery Navier-Stokes code Glenn-HT is coupled 
to a 3-D BEM steady state heat conduction solver. Glenn-HT is a multi-block cell-centered finite 
volume explicit code using a multi-stage Runge-Kutta based multigrid method time marching. 
The steady-state solution is sought by marching in time until dependent variables reach their 
steady-state values. The steady heat conduction equation is solved using the BEM with 
isoparametric bilinear discontinuous elements. We choose to employ discontinuous elements as 
they provide high levels of accuracy in computed heat flux values without resorting to special 
treatment of comer points required by continuous elements particularly when first kind boundary 
conditions are imposed to the conduction solver as is the case in the algorithm adopted in this 
paper. Moreover, the use of discontinuous elements throughout the BEM model eliminates much 
of the overhead associated with continuous elements, in particular, there is no need to generate, 
store, or access a connectivity matrix when using discontinuous elements. Details of the 
interpolation used to exchange nodal temperature and flux information from the disparate CFD 
and BEM grids are discussed. Results from a CHT numerical simulation of a 3-D film-cooled 
blade section are presented and are compared with those obtained from the standard approach of 
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a two temperature model. A significant difference in the level and distribution of the metal 
temperatures is found between the two models. 

Finally, current developments of an iterative strategy accommodating large numbers of 
unknowns by an artificial subsectioning of the blade are presented. The blade is subsectioned in 
the spanwise direction and a specially tailored iterative scheme is developed to solve the 
conduction problem with each subsection BEM problem solved using a direct LU solver. An 
adiabatic intial guess may be made for the sub-structure interface BEM nodes. Although the 
iterative method converges in some cases the iteration may be slow to converge. A better initial 
guess is provided via a physically-based initialization of the substructure interfacial temperatures. 
This is shown to provide an effective starting point for the iterative algorithm and significantly 
reduce the number of iterations required to achieve convergence. Results from several 
simulations in 2-D and 3-D show the process converges efficiently and offers substantial 
computational and storage savings. 
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1 Introduction 


Engineering analysis of complex mechanical devices such as turbomachines requires 
ever-increasing fidelity in numerical models upon which designers rely in their efforts to attain 
demanding specifications placed on efficiency and durability of modem machinery. 
Consequently, the trend in computational mechanics is to adopt coupled- fie Id analysis to obtain 
computational models which attempt to better mimic the physics under consideration, see Kassab 
and Aliabadi [1], The coupled field problem which we address in this paper is conjugate heat 
transfer (CHT): the coupling of convective heat transfer external to the solid body of a thermal 
component coupled to conduction heat transfer within the solid body of that component, see 
Figure 1. Conjugate heat transfer thus applies to any thermal system in which multi-mode 
convective/conduction heat transfer is of particular importance to thermal design, and thus CHT 
arises naturally in most instances where external and internal temperature fields are coupled. 



Figure 1. CHT problem: external convective heat transfer coupled to heat conduction within the 
solid. 


Conjugacy is often ignored in most analytical solutions and numerical simulations. For 
instance, it is common practice in analysis of turbomachinery, Heidmann et al. [2], to carry out 
separate flow and heat conduction analyses. Heat transfer coefficient as well as film effectiveness 
values are predicted using two independent external flow solutions each computed by imposing a 
different constant wall temperature at the surfaces of the turbine blade exposed to hot gases and 
film cooling air. The film effectiveness determines the reference temperature for the computed 
film coefficients. In turn, these values are used to impose convective boundary conditions to a 
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conduction solver to obtain predicted metal temperatures. As will be shown in the example 
section of this paper, the shortcomings of this approach which neglects the effects of the wall 
temperature distribution on the development of the thermal boundary layer are readily overcome 
by a CHT analysis in which the coupled nature of the field problem is explicitly taken into 
account in the analysis. 

There are two basic approaches to solving coupled field problems. In the first approach, a 
direct coupling is implemented in which different fields are solved simultaneously in one large 
set of equations. Direct coupling is mostly applicable for problems where time accuracy is 
critical, for instance, in aero-elasticity applications where the time scale of the fluid motion is on 
the same order as the structural modal frequency. However, this approach suffers from a major 
disadvantage due to the mismatch in the structure of the coefficient matrices arising from 
BEM, FEM and/or FVM solvers. That is, given the fully populated nature of the BEM 
coefficient matrix, the direct coupling approach would severely degrade numerical efficiency of 
the solution by directly incorporating the fully populated BEM equations into the sparsely banded 
FEM or FVM equations. A second approach which may be followed is a loose coupling strategy 
where each set of field equations is solved separately to produce boundary conditions for the 
other. The equations are solved in turn until an iterated convergence criterion, namely continuity 
of temperature and heat flux, is met at the fluid-solid interface. The loose coupling strategy is 
particularly attractive when coupling auxiliary field equations to computational fluid dynamics 
codes as the structure of neither solver interferes in the solution process. 

Several approaches can be taken to solve coupled field problems and most are based on 
either finite elements (FEM) or finite volume methods (FVM), or a combination of these two 
field solvers. Examples of such loosely coupled approaches applied to a variety of CHT problems 
ranging from engine block models to turbomachinery can be found in Comini et al. [3], Shyy and 
Burke [4], Patankar [5], Kao and Liou [6], Hahn et al. [7], Bohn et al. [8,9], and in Tayla et al. 
[10] where multi-disciplinary optimization is considered for CHT modeled turbine airfoil 
designs. Hassan et al. [11] develop a conjugate algorithm which loosely couples a FVM-based 
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Hypersonic CFD code to an FEM heat conduction solver in an effort to predict ablation profiles 
in hypersonic re-entry vehicles. Here, the structured grid of the flow solver is interfaced with the 
un-structured grid of heat conduction solvers in a quasi-transient CHT solution tracing the re- 
entry vehicle trajectory. Issues in loosely coupled analysis of the elastic response of solid 
structures perturbed by external flowfields arising in aero-elastic problems can be found in 
Brown [12] and Dowell and Hall [13]. In either case, the coupled field solution requires complete 
meshing of both fluid and solid regions while enforcing solid/fluid interface continuity of fluxes 
and temperatures, in the case of CHT analysis, or displacement and traction, in the case of aero- 
elasticity analysis. 

A different approach taken by Li and Kassab [14,15] and Ye et al. [16] who develop a 
BEM-based CHT algorithm thereby avoiding meshing of the solid region for the conduction 
solution. The method couples the boundary element method (BEM) to a FVM Navier-Stokes 
solver and was applied to solving two-dimensional steady state compressible subsonic CHT 
problems over cooled and uncooled turbine blades. The conduction problem requires solution of 
the Laplace equation for the temperature (or the Kirchhoff transform in the case of temperature 
dependent conductivity), and, as such, only requires a boundary discretization thereby 
eliminating the onerous task of grid generation within intricate regions of the solid. The boundary 
discretization utilized to generate the computational grid for the external flow-field can be 
considerably coarsened to provide the boundary discretization required for the boundary element 
method. Most modern grid generators used in computational fluid dynamics, for instance, 
GridPro™ [17], the topology-based algebraic grid generator used in the examples presented in 
this paper, allow for the multigrid option. Several levels of coarse discretization can thus be 
readily obtained. Furthermore, the BEM/FVM method offers the additional advantage of 
providing heat flux values and this stems from the fact that nodal unknowns which appear in the 
BEM are the surface temperatures and heat fluxes. Consequently, solid/fluid interfacial heat 
fluxes which are required to enforce continuity in CHT problems are naturally provided by the 
BEM conduction analysis. This is in sharp contrast to domain meshing methods such as FVM 
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and FEM where heat fluxes are computed by numerical differentiation in a post-processing stage. 
He et al. [18,19] adopted the BEM/FVM approach in further studies of CHT in incompressible 
flow in ducts subjected to constant wall temperature and constant heat flux boundary conditions. 
Kontinos [20] also adopted the BEM/FVM coupling algorithm to solve the CHT over metallic 
thermal protection panels at the leading edge of the X-33 in a Mach 15 hypersonic flow regime. 
Rahaim and Kassab [21] and Rahaim et al. [22] adopt a BEM/FVM strategy to solve time- 
accurate CHT problems for supersonic compressible flow over a 2-D wedged, and they present 
experimental validation of this CHT solver. In their studies, the dual reciprocity BEM [23] was 
used for transient heat conduction, while a cell-centered FVM was chosen to resolve the 
compressible turbulent Navier-Stokes equations. 

In this paper, we report on progress in the development and application of a BEM-based 
temperature forward/flux back (TFFB) coupling algorithm developed to solve conjugate heat 
transfer (CHT) arising in 3-D film-cooled turbine blades. The NASA-Glenn turbomachinery 
Navier-Stokes code Glenn-HT is coupled to a 3-D BEM steady state heat conduction 
solver[24,25]. The steady-state solution is sought by marching in time until dependent variables 
reach their steady-state values, and, as such, intermediate temporal solutions are not physically 
meaningful. In this mode of solving the steady-state problem, time-marching can be viewed as a 
relaxation scheme, and local time-stepping and implicit residual smoothing are used to accelerate 
convergence. The steady heat conduction equation reduces to the Laplace equation, and it is 
solved using the BEM with isoparametric bilinear discontinuous elements. We choose to employ 
discontinuous elements as they provide high levels of accuracy in computed heat flux values 
especially at sharp corners regions where first kind boundary conditions are imposed without 
resorting to special treatment of corner points required by continuous elements in particular when 
first kind boundary conditions are imposed [26,27]. In this application, sharp corners occur in 
many locations and first kind boundary conditions are imposed on all metal surfaces. Moreover, 
the use of discontinuous elements throughout the BEM model eliminates much of the overhead 
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associated with continuous elements, in particular, there is no need to generate, store, or access a 
connectivity matrix when using discontinuous elements. 

In order to resolve the flow physics, the CFD grid must be clustered in many regions. The 
BEM grid does not require such fine clustering and consequently the two grids are of quite 
different coarseness. The details of the interpolation used to exchange nodal temperature and flux 
information from the disparate CFD and BEM grids are presented. Results from a CHT 
numerical simulation of a 3-D film-cooled blade section are presented and results are compared 
with those obtained from the standard approach of a two temperature model. Significant 
difference in the level and distribution of the metal temperature is found between the two- 
temperature and CHT models. Finally, in order to address the large number of unknowns 
appearing in the 3-D BEM model, current developments of a strategy of artificial subsectioning 
of the blade are presented. Here, the approach is to subsection the blade in the spanwise 
direction. A specially tailored iterative scheme is developed to solve the conduction problem 
with each subsection BEM problem solved using a direct LU solver. A physically based initial 
guess is used to provide a good starting point for the iterative algorithm. Results from 2-D and 
3-D simulations show the process converges efficiently and offers substantial computational and 
storage savings. 

2 Governing Equations 

We first present the governing equations for the coupled field problem under 
consideration. The CHT problems arising in turbomachinery involves external flow fields that 
are generally compressible and turbulent, and these are governed by the compressible Navier- 
Stokes equations supplemented by a turbulence model. Heat transfer within the blade is governed 
by the heat conduction equation. Linear as well as non-linear options are considered. However, 
fluid flows within internal structures to the blade, such as film cooling holes and channels, are 
usually low-speed and incompressible. Consequently, density-based compressible codes tend to 
experience numerical difficulties in modeling such flows, unless low Mach number pre- 
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conditioning is implemented, see Turkel [28,29]. The Glenn-HT code is specialized to 
turbomachinery applications for which air is the working fluid and which is modeled as an ideal 
gas. 


2.1 Governing Equations for The Flow Field 

The governing equations for the flow field are the compressible Navier-Stokes equations, 
which describe the conservation of mass, momentum and energy. These can be written in integral 
form as 

F - j) • ndT = J Sdn (1) 

nr n 


'dW 

~dt 


dto + 


where H denotes the volume, T denotes the surface bounded by the volume O, and n is the 
outward-drawn normal. The conserved variables are contained in the vector W = (p. pu, pv, 
pw, pe, pk, puS) , where, p, u, v, w, e, k,u are the density, the velocity components in x-, y-, 
and z-di rections, and the specific total energy. The kinetic energy of turbulent fluctuations is 
denoted by k and the specific dissipation rate c o appear in the two equation k-u Wilcox 
turbulence model [30, 31] with modifications by Menter [32] and Chima [33] as implemented in 
Glenn-HT. The vectors F and T are convective and diffusive fluxes respectively, S is a vector 
containing all terms arising from the use of a non-inertial reference frame as well as production 
and dissipation of turbulent quantities. The working fluid is air, and it is modeled as an ideal gas. 
A rotating frame of reference can be adopted for modeling of rotating flows. The effective 
viscosity is given by 

P = Pi + Pt (2) 


where = pk/ui. The thermal conductivity of the fluid is then computed by a Prandtl number 
analogy where 


k f 


7 Pi Pt 

7 — 1 [ Pti Prt 


(3) 
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and Pr is the Prandtl number and 7 is the specific heat ratio. The subscripts l and t refers to 
laminar and turbulent values respectively. 

2.2 The Governing Equations of the Heat Conduction Field 

In the steady-state CHT solutions obtained in this paper, the NS equations are solved to 
steady state by a time marching scheme converging towards steady-state. A steady heat 
conduction analysis is carried out using the BEM at each time level chosen for the external flow- 
field and internal conduction field to interact in the iterative process. As such, the governing 
equation under consideration is 

V • [ k(T s )VT s } = 0 (4) 

where, T s denotes the temperature of the solid, and k s is the thermal conductivity of the solid 
material. If the thermal conductivity is taken as constant, then the above reduces to the Laplace 
equation for the temperature. When the thermal conductivity variation with temperature is an 
important concern, the nonlinearity in the steady-state heat conduction equation can readily be 
removed by introducing the classical Kirchhoff transform, U (T), see [34-37], which is defined 
as 

1 f T 

U (T)=— / k s (T)dT (5) 

k 0 Jt 0 

where T 0 is the reference temperature and k 0 is the reference thermal conductivity. The 
transform and its inverse are readily evaluated, either analytically or numerically, and the heat 
conduction equation transforms to a Laplace equation for the transform parameter U(T). The 
heat conduction equation thus reduces to the Laplace equation in any case, and this equation is 
readily solved by the BEM. 

In the conjugate problem, continuity of temperature and heat flux at the blade surface, T, 
must be satisfied: 
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(6) 


Tf = T s 

, dTf _ 8T S 
kf dn ~ h dn 

Here, Tf is the temperature computed from the N-S solution, T s is the temperature within the 
solid which is computed from the BEM solution, and d/dn denotes the normal derivative. Both 
first kind and second kind boundary conditions transform linearly in the case of temperature 
dependent conductivity. In such a case, the fluid temperature is used to evaluate the Kirchhoff 
transform and this is used a boundary condition of the first kind for the BEM conduction solution 
in the solid. Subsequently the computed heat flux, in terms of U, is scaled to provide the heat 
flux which is in turn used as an input boundary condition for the flow-field. 

3 Field Solver Solution Algorithms 

A brief description of the Glenn-HT code is given in this section. Details of the code and 
its verification in turbomachinery application can be found in [2, 38-41]. The heat conduction 
equation is solved using BEM. 


3.1 Navier-Stokes Solver 

Glenn-HT uses a cell-centered FVM to discretize the NS equations. Equation (1), is integrated 
over a hexahedral computational cell with the nodal unknowns located at the cell center (z, j. k ). 
The convective flux vector is discretized by a central difference supplemented by artificial 
dissipation as described in Jameson et al. [42]. The artificial dissipation is a blend of first and 
third order differences with the third order term active everywhere except at shocks and locations 
of strong pressure gradients. The viscous terms are evaluated using central differences. The 
overall accuracy of the code is second order, Heidmann et al. [2]. The resulting finite volume 
equations can be written at every computational node as 


dW 


Vi 


i,j,k ' 


i,j,k 


dt 


+ q 


i,j,k 


■ d ... = s ... 

~ i.j.k ~ i.j.k 


(7) 
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where W . . , is the cell-volume averaged vector of conserved variables, q and d ... are the 

~ & ~ i,j,k 

net flux and dissipation for the finite volume obtained by surface integration of Eqn. (1), and 
s . . , is the net finite source term. The above is solved using a time marching scheme based on a 
fourth order explicit Runge-Kutta time stepping algorithm. The steady- state solution is sought by 
marching in time until the dependent variables reach their steady-state values, and, as such, 
intermediate temporal solutions are not physically meaningful. In this mode of solving the 
steady-state problem, time-marching can be viewed as a relaxation scheme, and local time- 
stepping and implicit residual smoothing are used to accelerate convergence. A multigrid option 
is available in the code. The code also adopts a multi-block strategy to model complex 
geometries associated with film-cooled blade problems. Here, locally structured grids blocks are 
generated into a globally unstructured assembly. 

Glenn-HT adopts a k-u> turbulence model which integrates to the wall and does not 
require maintaining a specified distance from the wall as no wall functions are used. The 
computational grid is sufficiently fine near the wall to yield a ;y + value of less than 1.0 at the first 
grid point away from the wall. A constant value of 0.9 is taken for the turbulent Prandlt number 
in all heat transfer computations, while a constant value of 0.72 is used for the laminar Prandtl 
number. Moreover, the temperature variation of the laminar viscosity is taken as a 0.7 power law, 
see Schlichting [43], and c p is taken as constant. 

3.2 Heat Conduction Boundary Element Solution 

The heat conduction equation reduces to the same governing Laplace equation in the 
temperature or the Kirchhoff transform. In the boundary element method this governing partial 
differential equation is converted into a boundary integral equation (BIE), see [44-46], as 

C(0 T(0+ i T(x)q*(x, OdS(x)= <f q(x)T*(x, £)dS(x) ( 8 ) 

Js Js 
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where S(x ) is the surface bounding the domain of interest, £ is the source point, x is the field 
point, q(x) = — kdT/dn is the heat flux, T*(x, £) is the so-called fundamental solution, and 
q*(x,£) is its normal derivative with d/dn denoting the normal derivative with respect to the 
outward-drawn normal. The fundamental solution (or Green free space solution) is the response 
of the adjoint governing differential operator at any field point x due to a perturbation of a Dirac 
delta function acting at the source point £. In our case, since the steady state heat conduction 
equation is self-adjoint, we have 

kV 2 T*(x,0 = ~8(x,0 (9) 

Solution to this equation can be found by several means, see for instance Liggett and Liu [47], 
Morse and Feschbach [48], and Kellog [49], as 

T %x,0= -^ l n r (x,0 in 2-D (10) 

= , , — -r in 3-D 

Auk r{x, £) 

where r(x, £) is the Euclidean distance from the source point £. The free term C(£) can be shown 
analytically to be 

C(i)=i - k d Al^L dS(x) (11) 

■1 Sir) on 

Moreover, introducing the definition of the fundamental solution in the above, it can be readily 
be determined that C(£) is the internal angle (in degrees in 2-D and in steradians in 3-D) 
subtended at source point divided by 27rin 2-D and by An in 3-D when the source point £is on 
the boundary and takes on a value of one when the source point £ is at the interior. 
Consequently, the free term takes on values 1 > C(£) > 0. In the standard boundary element 
method, the BIE is discretized using two levels of discretization: 

1. the surface S is discretized into a series of j = 1, 2. ..TV elements. This is traditionally 
accomplished using polynomial interpolation, bilinear and biquadratic being the most 
common. In general, 
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( 12 ) 


N 

s = J>s, 

‘ 1=1 

and on each surface element A Sj the geometry is discretized using local shape 
functions iV fc (p, (') in terms a homogeneous coordinates (//. (') which each take on 
values between [ — 1, 1] as 

NGE 

Xjiv, 0 = (13) 

fc= 1 
NGE 

k= 1 
NGE 

zMO = 'E Nl( - r >.<)^ 

k = 1 

Here, Zj) denote the location of the k = 1, 2. ..NGE boundary nodes used to 

define the boundary element j geometry. 


2. the distribution of the temperature and heat flux is modeled on the surface. This is 
usually accomplished using polynomial interpolation as well. Common discretizations 
include: constant (where the mean value of T and q are taken on an element surface), 
bilinear, or biquadratic. In general, 

NPE 

T,G,0 = Y. Mk ^G)T^ (14) 

k = 1 

NPE 

q,(v, C) = £>*(,, C)?y 

k = 1 

It is noted that the order of discretization of the temperature and heat flux need not be 
the same as that used for the geometry, leading to subparametric (lower order than that 
used for the geometry), isoparametric (same order than that used for the geometry), 
and superparametric (higher order than that used for the geometry) discretizations. 
Moreover, the temperature and heat flux are discretized using k = 1,2, ...NPE 
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discrete nodal values whose location within the element j can be chosen to 

(a) coincide with the location of the geometric nodes: continuous elements. 

(b) be located offset from the geometric nodes: discontinuous elements. 

We choose to employ bilinear discontinuous isoparametric elements as they provide high 
levels of accuracy in computed heat flux values especially at sharp comers regions where first 
kind boundary conditions are imposed without resorting to special treatment of comer points 
required by continuous elements [26,27]. In this type of boundary element, the field variables T 
and q are modeled with discontinuous bilinear shape functions across each element while the 
geometry is represented locally as continuous bilinear surfaces. Figure 2 below shows a typical 
bilinear isoparametric boundary element along with its transformed representation in the local 
coordinate rj-( system. 



Figure 2. Bilinear isoparameteric discontinuous boundary element. 


The global coordinate system (x, y, z) is transformed into a local coordinate system ( 77 , () by 
the backward transformation 
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( 15 ) 


Xj(v, 0 = ^N k {r],C)x k j 

k = 1 

yj(v, 0 = 

fc=i 

z j (ri,0='£ Nk ^-0^ 

k = 1 

where (x k ,y k ,Zj) are the global positions of the /c-th geometric node of element j. The four 
bilinear shape functions are defined as 

N \v, 0 = ^(i -v)0- -C) 

N 2 (v,C) = t(1 + »7)(1-C) 

f (16) 

7V 3 (? 7 ,C) = i (l + r?)(l + C) 
iV 4 (r 7 ,C) = ^(l-r 7 )(l + C) 


The field variables, T and q, are modeled to vary bilinearly across the boundary element through 
the use of four discontinuous shape functions with nodes located at an off-set position of 12.5% 
from the edges of the element. The field variables and shape functions are described as follows: 

4 

T ] (vX) = J2 Mi (’l,0T? (17) 

k = 1 
4 

®(>?,C ) = Y, Mt ^,Qq t i 

k = 1 


with the bilinear shape functions M k defined as, 

0 = ^(3 -477) (3 -4C) (18) 

do 

M 2 (v,C) = ^(3 + 4t/)(3-4C) 

do 

^ 3 (>l.C) = ^(3 + 4r f )(3 + 4C) 

do 

M 4 (r),() = ---(3 — 477 ) (3 + 4£) 
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Introducing the above discretization in the BIE in Eqn. (8), and collocating the discretized BIE at 
each of the boundary nodes there results 


N NPE N NPE 

Cte) m,) + E E = E E g «?‘ 

j=i fc=i i =1 fc =1 


( 19 ) 


where the influence coefficients if and G are defined as 

l J 


TT K 

— 


Gr — 


AS,- 


M k (r]X) dS(x) 

j 

T*(x, £i) M k (r], () dS(x) 


( 20 ) 


A5, 


These influence coefficients are evaluated numerically via Gauss-Legendre quadratures with 
special adaption when evaluating the integrals on AS) (for the element upon which the source 
point lies) and heuristic adaptive quadratures for elements that are close to the node of interest 
(see Appendix). The surface integrals in the above equation depend purely on the local geometry 
of the element and the location of the source point Upon collocation of the above at every 
boundary node where the temperature and heat flux are defined, the following algebraic form is 
obtained: 

[H\ {T,} = [G]{<&} (21) 


Here the influence matrices [H] and [G] are evaluated numerically using quadratures. Once the 
boundary conditions are specified, the above is re-arranged in standard [A] {a:} = {b} and the 
ensuing equations are solved by direct or iterative methods. In a fully conjugate solution using 
the algorithm described in this paper, these BEM equations are solved subject to the following 
boundary condition at external and internal bounding walls which are in contact with the fluid 
and denoted by T conjugate: 


T s 


sir 


conjugate 


T f 


( 22 ) 
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In the reduced periodic 3-D computational model to be discussed in the example section, 
adiabatic conditions are also imposed at the flowfield periodic surfaces in the spanwise direction, 
i.e. there 

q s = 0 (23) 

Once these equations are solved, the heat flux is known at all surface nodes. This is the sought- 
after quantity in the CHT algorithm to be shortly outlined. In the case where the conduction 
problem is solved without further treatment, the basic BEM code had options of using an LU 
decomposition for small numbers of equations and a GMRES iterative solver with an incomplete 
LU (ILU) pre-conditioning for large numbers of equations. When the number of equations gets 
very large, storage becomes an important issue, as the coefficient matrix is fully-populated. We 
will discuss an effective treatment of such problems in a later section. 

3.3 CHT Algorithm 

The Navier-Stokes equations for the external fluid flow and the heat conduction equation 
for heat conduction within the solid are interactively solved to steady state through a time- 
marching algorithm. The surface temperature obtained from the solution of the Navier-Stokes 
equations is used as the boundary condition of the boundary element method for the calculation 
of heat flux through the solid surface. This heat flux is in turn used as a boundary condition for 
the Navier-Stokes equations in the next time step. This procedure is repeated until a steady-state 
solution is obtained. In practice, the BEM is solved every few cycles of the FVM to update the 
boundary conditions, as intermediate solutions are not physical in this scheme. In the calculations 
carried out in this study, the BEM solution was run every ten cycles of the finite volume solver. 
This is referred to as the temperature forward/flux back (TFFB) coupling algorithm as outlined 
below: 
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FVM Navier-Stokes solver: 


1. Begins with initial adiabatic boundary condition at solid surface. 

2. Solves compressible NS for fluid region. 

3. Provides temperature distribution to BEM conduction solver 
after a number of iterations. 

4. Receives flux boundary condition from the BEM as input for 
next set of iterations. 

• BEM conduction solver: 

1. Receives temperature distribution from FVM solver. 

2. Solves steady-state conduction problem. 

3. Provides flux distribution to FVM solver. 

The transfer of heat flux from the BEM to the FVM solver is accomplished after under- 
relaxation. 

<1 = P - P)Cl M (24) 


with /3 taken as 0.2 in all reported calculations. The choice of the relaxation parameter is through 
trial and error. In certain cases, it has been our experience that a choice of larger relaxation 
parameter can lead to non-convergent solutions [50]. The process is continued until the NS solver 
converges and wall temperatures and heat fluxes converge, that is until Eqn. (6) is satisfied 
within a set tolerance 


T .-T 

r^j J r^j S 


< ex 


Q 


f 


- q 


< 6q 


(25) 


where the tolerances ex and e q are taken as 0.001. 

It should be noted that alternatively, the flux could be specified as a boundary condition 
for the BEM code leading to a flux forward temperature back (FFTB) approach. However, when 
a fully conjugate solution is undertaken, this would amount to specifying second kind boundary 
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conditions completely around the surface of a domain governed by an elliptic equation, resulting 
in a non-unique solution. The TFFB algorithm avoids such a situation. 

3.4 Interpolation Between BEM and FVM Grids 

An issue arises in information transfer between CFD and BEM as there exists a 
significant difference in the levels of discretization between the two meshes in a typical CHT 
simulation. Accurate resolution of the boundary layer requires a FVM surface grid which is much 
too fine to be used directly in the BEM. A much coarser surface grid is typically generated for the 
BEM solution of the conduction problem. The disparity between the two grids requires a general 
interpolation of the surface temperature and heat flux between the two solvers as it is not 
possible in general to isolate a single BEM nodes and identify a set of nearest FVM nodes. 
Indeed in certain regions where the CFD mesh is very fine, a BEM node can readily be 
surrounded by tens or more FVM nodes. 

A distance-weighted interpolation, reminiscent of radial basis function (RBF) 
interpolation [23], is adopted for transfer of temperature and flux values between the BEM and 
CFD grids. Consider Fig. 3(a), here the location of a BEM node is identified on the right-hand- 
side by a star-like symbol. Let us consider the problem of transferring the temperature from the 
FVM grid to the BEM grid. Let us denote the position of the BEM node of interest by f u and the 
location of an FVM node by r :r The radial distance from every FVM node to the BEM node of 
interest is then r ?y = |r, — r t \. Let us suppose that the number of all FVM surface nodes lying 
within a ball of radius R ma x centered 

about r is N} )a u. Moreover, let us denote two cases. Case I where all r tJ > e and case II where 
there is an FVM node located at r, if such that r, ? < e, where e is a tolerance. Then, the value of 
the temperature at the BEM node r 3 is evaluated as 
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TBEhlifi) 


for case I 


(26) 
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for case II 


In all calculations, the maximum radius R max of the sphere is set to 2.5% of the maximum 
distance within the solid region and e is set to e = R max x 10 20 ■ These limits may be adjusted to 
suit the problems at hand. 


T 


q 

(a) transfer of nodal temperatures and fluxes between CFD and BEM grids. 



(b) 2-D illustration of five CFD nodes nearest to a BEM node located at r, . 

Figure 3. Transfer of nodal values from FVM and BEM (and back) independent surface meshes 
is performed with a distance weighted radial interpolation. 
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4 Strategy for BEM Models of Large-Scale Three-Dimensional Heat Conduction Problems 

As mentioned, the BEM is ideally suited for the solution of linear and non-linear heat 
conduction problems and is a particularly advantageous numerical method due to its boundary- 
only feature, however, the coefficient matrix of the resulting system of algebraic equations is 
fully populated. For large-scale problems that occur in engineering 3-D modeling of complex 
structures this poses very serious numerical challenges due to large storage requirements and 
iterative solution of large sets of non-sparse equations. This problem has been approached in the 
BEM community by one of two approaches: the artificial subsectioning of the 3-D model into a 
multi-region model in conjunction with block-solvers reminiscent of the FEM frontal solvers 
[51,52] and the adoption of multipole methods in conjunction with the GMRES non-symmetric 
iterative solver [53-55]. The first approach of domain decomposition or subsectioning produces 
a sparse block coefficient matrix that is efficiently stored and has been successfully implemented 
in commercial codes such as BETTI and GPBEST in the context of continuous boundary 
elements. However, the method requires generation of complex data- structures identifying 
connecting regions and interfaces prior to analysis. The second approach is very efficient, 
however, it requires complete re-write of the BEM code to adopt multipole formulation. 
Recently, a novel technique using wavelet decomposition has been recently proposed to reduce 
matrix storage requirements without need for a major alteration of traditional BEM codes [56]. 

We propose to adopt the first approach, however, we do not use a block solver but rather 
a region-by-region iterative solver. Although it was reported earlier in the literature that this 
process sometimes has difficulty converging for nonlinear problems [34,35], it has been shown 
[36] that when properly implemented, the iterative process converges very efficiently and can 
offer very substantial savings in memory. Moreover, the technique does not require any complex 
data-structure preparation. Indeed, the approach is somewhat transparent to the user, a significant 
advantage in coupling BEM to other field solvers. It should be noted that this sub- sectioning 
method is under current development and has not yet been integrated into the CHT solver at the 
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point of writing this paper, and thus the technique along with an example 3-D conduction 
solution are presented herein with this explicit caveat. 

The solution algorithm for the multi-region BEM iteration process consists on the steps 
described as follows. First, the problem domain Q is identified along with the corresponding 
boundary conditions over the boundary T. A typical problem definition along with the 
corresponding boundary conditions and a sample single-region BEM discretization is depicted 
below. 


q=h 4 (T-T ao4 ) r 4 



q-h 2 (T T^t) r 2 


Figure 4. BEM problem domain, boundary conditions, and single region discretization. 

If a standard BEM solution process were to be adopted a system of influence coefficient 
matrices and boundary values of size N, where N is the number of boundary nodes used to 
discretize the problem, will be formed. The number of floating point operations required to arrive 
at the algebraic system is proportional to N 2 as well as direct memory allocation also 
proportional to N 2 . With the aid of the boundary condition distribution, the system is re- 
arranged as 

mm = [G]M =c [A]{x} = { b } (27) 

where {rc} represents the unknowns {T} or {q} around the boundary nodal distribution. The 
solution to the algebraic system for the boundary unknowns can be performed using a direct 
solution method such as LU decomposition requiring floating point operations proportional to 
N 3 or an indirect method such as Bi-conjugate Gradient or General Minimization of Residuals 
which, in general, require floating point operations proportional to N 2 to achieve convergence. 
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If a multi-region BEM iteration solution process were to be adopted instead, the domain 
is divided into K subdomains and each one is independently discretized. It is worth mentioning 
that the BEM discretizations of neighboring subdomains does not have to be coincident, this is, 
at the connecting interface, boundary elements and nodes from the two colliding subdomains are 
not required to be structured following a sequence or particular position, the only requirement at 
the connecting interface is that it forms a closed boundary with the same path on both sides. 
Later, it will be shown that the information between neighboring subdomains separated by an 
interface will be passed through an interpolation process as opposed than just a node to node 
connection. 


q=h 4 (T- Tj 
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Figure 5. BEM problem domain, boundary conditions, and single region discretization. 


Figure 5 above shows the same problem depicted in Fig. 4 with a multi-region BEM 
discretization of four ( K = 4) subdomains. The boundary value problem will now be solved 
independently over each subdomain where initially, a guessed boundary condition is imposed 
over the interfaces in order to ensure the well-posedness of each subproblem. For instance, the 
boundary value problem of subdomain is transformed into the algebraic analogue of 
corresponding influence coefficient matrices and nodal boundary values as 

W 2 T ni (x,y) = 0 =>- [H Q] ] {Tr, } = (28) 

The composition of this algebraic system requires floating point operations proportional to the 
square of the number boundary nodes in the subdomain (n 2 )as well as for direct memory 
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allocation (n 2 ). This new proportionality number n is roughly equivalent to 


n 


2N 
K + l 


( 29 ) 


as long as the discretization along the interfaces has the same level of resolution as the 
discretization along the boundaries. Direct memory allocation requirement for later algebraic 
manipulation is now reduced to a proportion of n 2 as the influence coefficient matrices can easily 
be stored in ROM memory for later use after the boundary value problems on remaining 
subdomains have been effectively solved. For the example shown here where the number of 
subdomains is K = 4 the new proportionality value n is approximately equal to n ~ 2N / 5. This 
simple multi-region example reduces the memory requirements to about 
n 2 / N 2 = (4/25) = 16% of the standard BEM approach. 

The algebraic system for subdomain Q\ is rearranged, with the aid of given and guessed 
boundary conditions, as 

[Hn.KToJ = \Co.]{qn } =* [A a j{x ni } = {boA (30) 

The solution of the new algebraic system of subdomain Q\ requires now a number of floating 
point operations proportional to n 3 / N 3 = (8/125) = 6.4% of the standard BEM approach if a 
direct algebraic solution method is employed, or a number of floating point operations 
proportional to n 2 / N 2 = (4/25) = 16% of the standard BEM approach if an indirect algebraic 
solution method is employed. For both, floating point operation count and direct memory 
requirement the reduction is dramatic. However, as the first set of solutions for the subdomains 
were obtained using guessed boundary conditions along the interfaces, the global solution needs 
to follow an iteration process and convergence criteria. 

Globally, the floating point count for the formation of the algebraic setup for all K 
subdomains must be multiplied by K , therefore, the total operation count for the coefficient 
matrices computation is given by, 
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4 K 


(31) 


K 


n 


N 2 


(K + l) 2 


or, for this particular case with K = 4, Kn 2 /N 2 = 16/25 = 64% of the standard BEM 
approach. Moreover, the more significant reduction is revealed in the RAM memory 
requirements as only the memory needs for one of the subdomains must be allocated at a time. 
The rest of the coefficient matrices for the remaining subdomains are temporarily stored in ROM 
memory until access and manipulation is required. Therefore, for this case of K = 4, the true 
memory reduction is n 2 / N 2 = 4/25 = 16% of the standard BEM approach. 

With respect to the algebraic solution of the system of Eqn. (30), if a direct approach as 
the LU decomposition method is employed for all subdomains, the LU factors of the coefficient 
matrices for all subdomains can be computed only once at the first iteration step and stored in 
ROM memory for later use during the iteration process for which only a forward and a backward 
substitution will be required to solve the system at hand. This feature allows a significant 
reduction in the operational count through the iteration process until convergence is achieved, as 
only a number of floating point operations proportional to n as opposed to n 3 is required at each 
iteration step. To this computation time is added the access to ROM memory at each iteration 
step which is usually larger than the access to RAM. Moreover, when applied to the CHT 
problem, again storage of the LU factors offers significant computational savings over using an 
iterative method. However, use of direct solvers requires each sub-section to be kept at a 
discretization close to 1,000 bi-linear boundary elements. 

The iteration process follows the initial step of guessing the interface conditions. This is a 
crucial step as the more physical information the initial guess incorporates the less iterations will 
be required to reduce the error. The simplest choice is to assume adiabatic conditions at the 
artificial interfaces. Results from several numerical studies show this approach leads to intial 
temperature fields that are far from the final temperature field and which are slow to update 
iteratively. A rather more efficient initial guess can be made using a physically based 1-D heat 
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conduction argument for every node on the external surface to every node at the interface. The 
following algebraic initial guess for any interfacial node can be readily derived, see [57] for 
details, as 


Ti 


3 = 1 


Nt N Nh 

E BijTj - E BijRijtj + E NN 

3 = 1 


3 = 1 


Hi i+l 


n t N h 

Si-E^ + Eitft 

3 = 1 J =1 


(32) 


where TVy, JV g , and TV/, are the number of first, second, and third kind boundary conditions 
specified at the external (non-interfacial) surfaces and 


D ''b • "J 
H,j = — (rij ■ rtj) 

J* A. 

3=1 |r «! 


(33) 


with N = Nt + N q + Nh , the thermal conductivity of the medium is k, the film coefficient at 
the j-th convective surface is h :n the outward-drawn normal to any surface is n :n the position 
vector from the interfacial node i to the external surface node j is r tJ and its magnitude is 
r,j = \r t j\, while the area of element j denote by Aj is readily computed as 

Aj = (p dT(x, y, z) = / Jj(riX)drid( (34) 


Once the initial temperatures are imposed as boundary conditions at the interfaces, a 
resulting set of normal heat fluxes along the interfaces will computed. These are then non- 
symmetically averaged in an effort to match the heat flux from neighboring subdomains. 
Considering a two-domain substructure the averaging at the interface is explicitly given as, 



gn, + qL 2 
2 


(35) 
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and, 



Qa 2 + gfl! 
2 


(36) 


to ensure the flux continuity condition after averaging. Compactly supported radial 

basis interpolation can be employed for the flux average to account for unstructured grids along 
the interface from neighboring subdomains. 

Using these fluxes the BEM equations are again solved leading to mismatched 
temperatures along the interfaces for neighboring subdomains. These temperatures are 

interpolated, if necessary, from one side of the interface to the other side using a compactly 
supported radial basis functions to account for the possibility of interface mismatch between the 
adjoining substructure grids. Once this is accomplished, the temperature is averaged out at each 
interface. Illustrating this for a 2 domain substructure, again we have for regions 1 interface, 




and region 2 interface 




(37) 


(38) 


In case a real or physical interface exists and a thermal contact resistance is present between the 
connecting subdomains, the temperatures are averaged out as, 



T 1 

1 Q 1 


+ + dV 

7 ) r n q Ql 


(39) 


and, 




+ R 


(40) 


where R" is the value of the thermal contact resistance imposing a jump on the interface 
temperature values. These now matched temperatures along the interfaces are used as the next set 
of boundary conditions. 
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Figure 6. Iteration along the interface (a) imposed heat flux, (b) mismatched resulting 
temperatures, and (c) interpolated and averaged out temperatures as new boundary conditions. 

Finally, in order to provide an improved initial guess to the BEM solver, a two-level 
discreetization approach is employed. The initial guess provided by physically-based procedure 
in Eq. (32) is used to solve the conduction solution for a constant element conduction model. 
This coarse computation is carried out using the same discretization as the bi-linear model, 
however, all nodes are collapsed to a single central node, and this requres l/64th of the work 
required to solve the bi-linear model. Upon convergence of the constant element model, a full bi- 
linear solution is subsequnetly computed with the constant element model providing the initial 
guess. 

This iteration process is continued until a convergence criterion is satisfied. A measure of 
convergence may be defined as the T 2 norm of mismatched temperatures along all interfaces as: 



This norm measures the standard deviation of BEM computed interface temperatures T 1 and 
averaged-out updated interface temperatures T ? f. The iteration routine can be stopped once this 
standard deviation reaches a small fraction e of AT, where AT is the maximum temperature 
span of the global field. This concludes the numerical developments in this paper, attention is 
now given to numerical examples. 
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5 Numerical Results and Discussion 


We now present results of a full conjugate solution of a film-cooled blade under operating 
conditions which match a planned experiment at NASA Glenn Research center and assumes 
periodicity in the spanwise direction for one pitch of film-cooling hole patterns. We compare 
results of this simulation to those obtained from the standard two temperature method. This 
simulation uses the standard BEM approach to heat conduction. We also present results from 
several simulations fo the sub-sectioning BEM approach to heat conduction modeling in 2-D and 
in 3-D. In 3-D, a square cross-section bar, a thrust vector control vane, and a cooled turbine vane 
geometry are used to illustrate the sub-sectioning method described in this report. 

5.1 CHT simulation of a 3-D Film-cooled Turbine Blade 

Film cooling is commonly used in turbine designs to produce a buffer layer of relatively 
cool air between the turbine blade and the hot freestream gas in the first and second rows of 
blades and vanes. The CHT computation is carried out on a computational model of a realistic 
film-cooled turbine vane accounting for the three-dimensional vane geometry including plena 
and film holes and is based on a Honeywell film-cooled engine design, see Heidmann et al. [2]. 
The geometry of this test vane is based on the engine vane midspan coordinates, and is scaled up 
by a factor of 2.943 to allow matching of engine exit Mach number (0.876) and exit Reynolds 
number (2.9 x 10 6 based on true chord) with atmospheric inlet conditions. The test vane has a 
true chord of 0.206 m. Since the test vane is of constant cross section, only one spanwise pitch 
of the film hole pattern was discretized, with periodicity of the flow-field enforced at each end. 
This simplification assumes no effect of endwalls, but greatly reduces the number of grid points 
required to model the vane. However, the thermal boundary conditions enforced at these ends in 
the conduction analysis were adiabatic. The vane has two plena which feed 12 rows of film 
cooling holes as well as trailing-edge ejection slots, see Fig. 7. Trailing edge ejection is blocked 
in the computation as the planned experiment has no slot cooling. Detailed geometrical data for 
each row of film holes as well as hole distribution are provided in [2]. A multi-block grid 
approach is adopted to model this complex geometry and generated the FVM grid using the 
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topology-based algebraic grid- generation program GridPro™ [17] with the final grid consisting 
of 140 blocks and a total of 1.2 x 10 6 finite volume computational cells. The FVM grid consists 
of 20 cells across both the inlet and outlet boundaries, 60 cells on the periodic boundary, over 
200 cells around the vane, and 44 cells from the vane to the periodic boundary. 



Figure 7. Film-cooled blade profile used in the CHT simulation. 


A blade-to-blade view of the FVM grid is shown in Fig. 8 and Fig. 9 shows the FVM grid in the 
leading edge region of the vane. 
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Figure 8. Blade to blade computational grid cross section. 



Figure 9. FVM grid in the leading edge region of the blade. 
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The flow conditions for all simulations use a free-stream inlet flow to the vane at an angle 
of 0° to the axial direction, with all temperatures and pressures normalized by the inlet stagnation 
values of 31097? and 10 atmospheres, respectively. The inlet turbulence intensity is set at 8.0% 
and the turbulence scale is 15.0% of vane true chord. Other inflow quantities are set by means of 
the upstream-running Riemann invariant. The vane downstream exit flow is defined by imposing 
a constant normalized static pressure of 0.576, which was empirically determined to yield a 
desired exit Mach number of 0.876. Periodicity was enforced in both the blade-to-blade and 
spanwise directions based on vane and film hole pitches, respectively. Moreover, in order to 
maintain a true periodic solution, inflow to the plena was provided by defining a region of each 
plenum wall as an inlet and introducing uniform flow normal to the wall. In Fig. 8, these regions 
are shown to lie on either side of the internal wall that separates the two plena. In practice, there 
will be spanwise flow in the plenum, but bleed of the plenum flow into the film holes results in a 
spanwise-varying mass flow rate and static pressure, which would violate spanwise periodicity 
imposed in this particular reduced computational model. The non-dimensionalized inflow 
stagnation temperature to the plena was 0.5, corresponding to a coolant temperature of 1554.57?. 
the velocity was fixed to the constant value required to provide the design mass flow rate to each 
plenum, and static pressure was extrapolated from the interior. The inflow patch for each 
plenum was defined to be sufficiently large to yield very low inlet velocities (Mach 
number < 0.05), allowing each plenum to approximate an ideal plenum. All solid walls were 
imposed with a no-slip boundary condition. The blade metal material is taken as Inconel with a 
conductivity of kuadte = 1-34 Btu/hrin R taken at 2174.97? which is estimated to be the average 
blade temperature. 

The FVM metal surface grid consists of 38, 000 cells at the 4th level of multi-grid. The 
grid was coarsened to generate a BEM grid of 1,3000 bilinear cells with 5,2000 nodal 
unknowns. Two cases are computed in the numerical simulation in order to obtain the metal 
temperature: 
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1. The traditional two-temperature approach, whereby two different isothermal wall boundary 
condition extended to all wall surfaces, including the film hole surfaces and plenum surfaces. 
Two solutions were generated with constant wall temperatures T w of T ir .\ = 2174.9/? and 
T Wt 2 = 2485.6 R imposed on all blade surfaces. The flowfield was computed from the plena 
through the cooling holes and over the blade. The predicted wall heat fluxes at each node q” 
computed from each of these isothermal solutions were used to solve simultaneously for 
adiabatic wall temperature, T aw , and heat transfer coefficient, h, referenced to the computed 
adiabatic wall temperature, under the assumption that T aw and h are independent of the wall 
temperature. That is at each node we have 

Qw = h(Tw,l — Taw) /42 s ) 

C = h(T w , 2 - T aw ) 

In turn, these film coefficient and associated adiabatic wall distributions were used in the 
BEM to compute metal temperatures. 

2. A full CHT solution was carried out using the same grids and boundary conditions above 
except at the blade surface where conjugate conditions were imposed. The conjugate 
solutions converged in 1000 iterations with a BEM conduction calculation performed each 10 
FVM iterations. The BEM code was written as a subroutine to the Glenn-HT code and 
subroutines were coded to exchange information between the two codes in terms of the FVM 
and BEM grids as well as boundary condition information. The Glenn-HT code was modified 
to allow for non-isothermal boundary condition specification. 

All computations were performed at NASA Glenn Research Center on an SGI Origin 
2000 cluster with 32 processors. Flow computations were carried out and considered converged 
when residuals were driven below 10 5 . Results of the blade surface temperatures predicted by 
the simulations are shown in Fig. 10 for the CHT solution and in Fig. 11 for the two constant 
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temperature approach. The two temperature distributions are markedly different with a 
temperature span of AT = 1720 — 24207? across the surface of the blade while the CHT 
solution predicted a temperature span of AT = 1620 — 26207? across the blade. In addition to 
CHT computations predicting lower minimum (1007? colder) and higher maximum temperatures 
(2007? hotter), the distribution of cold and hot regions are quite different as is evident from the 
surface plots. For instance, with conduction taken into consideration in the CHT simulation, the 
thin trailing regions are seen to reach higher temperatures than predicted by the isothermal 
approach, while the forward plenum region is seen to be effectively cooler. This has severe 
implications in materials design and subsequent thermal stress analysis of the blade carried out 
using these metal temperatures. 



T: 1620 1720 1820 1920 2020 2120 2220 2320 2420 2520 2620 

Figure 10. Blade surface temperature predicted by the CHT solution. 
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Figure 1 1 . Blade surface temperature predicted by the BEM using h and T aw provided from the 
two temperature approach. 

5.2. Steady Conduction BEM Modeling Using Sub-Sectioning and Iteration 

Results are now presented for several simulations using the sub-sectioning iterative 
method for a pure heat conduction problem. In the first examples the initial guess is carried out 
using adiabatic conditions at the interface while in the final examples the physically-based 
reasoning is used to initiate the iteration. The parameter e is set to 10 3 . All computations were 
performed on a Pentium 4, 1.8 GHz PC with 512 MB 800MHz RDRAM. 

5.2.1 Steady Conduction Model of a 10-Region 2D Slab 

A 2D rectangular slab of dimension 10x1 is initially discretized in a single-region using 
440 equally spaced quadratic isoparametric discontinuous boundary elements with a total of 1320 
nodes around the boundary. Figure 12 shows the boundary conditions and the BEM discretization 
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for this problem. The isotherm distribution is plotted in Figure 13. The thermal conductivity was 
imposed as k = 1 W/mK. 



Figure 12. Boundary conditions and BEM discretization of single-region 2D rectangular slab. 



T: 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 



Figure 13. Isotherm distribution for single-region 2D rectangular slab. 

Next, the rectangular slab is divided into 10 sub-domains and discretized with 80 quadratic 
isoparametric discontinuous boundary elements with 240 nodes. Figure 14 shows the boundary 
conditions and the BEM discretization for this problem. The isotherm distribution is plotted in 
Figure 15 for the converged solution after 12 iterations with adiabatic conditions used as initial 
guess at the subsection interface nodes. 



Figure 14. Boundary conditions and BEM discretization of 10-region 2D rectangular slab. 
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T: 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 



Figure 15. Isotherm distribution for 10-region 2D rectangular slab. 

The plot of the L 2 norm progression is shown in Figure 16 for the first 20 iterations, however, 
the level of convergence of e • AT was achieved after 12 iterations. Table 1 shows the memory 
requirement proportions for each case and the computation time for the algebraic setup and 
solution after the 12 iterations for which L 2 = 0.04. 



Figure 16. L 2 norm progression for 10-region 2D rectangular slab. 


Table 1. Memory and time comparison for 2D rectangular slab problem. 



1 -Region 

10-Region 

Iterations 

0 

12 

L 2 Norm 

0 

0.04 

Memory 

100% 

3.3% 

Time for Setup 

50s 

16s 

Time for Solution 

64s 

50s 
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5.2.2 Steady Conduction model of a 5-Region 3D Rectangular Bar 

A 3D 5x1 rectangular bar is initially discretized in a single-region using 550 equally spaced 
bilinear isoparametric discontinuous boundary elements with a total of 2200 nodes around the 
boundary. Figure 17. a. shows the BEM discretization for this problem along with the isotherm 
distribution in Figure 17.b. The boundary conditions are distributed as first kind on both end 


faces with a temperature T = 0 6, adiabatic 
remaining three surfaces with = 100°C and 
imposed as k = 1 W /mK. 



on the bottom face, and convective on the 
h = 1 W /m 2 K. The thermal conductivity was 



T: 0 10 20 30 40 50 60 70 80 90 100 


(a) (b) 

Figure 17. Rectangular bar example: (a) single region BEM discretization and (b) Isotherm 
distribution of single-region 3D rectangular slab. 

Next, the 3D slab is divided into 5 sub-domains and discretized with 150 bilinear 
isoparametric discontinuous boundary elements with 600 nodes. Figure 18. a. shows the BEM 
discretization for this problem and the isotherm distribution is plotted in Figure 18.b. for the 
converged solution after 4 iterations. Adiabatic conditions are used as an initial guess at the sub- 
section interface nodes. 
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T: 0 10 20 30 40 50 60 70 80 90 100 


(a) (b) 

Figure 18. Five sub-section Model for a rectangular bar: (a) BEM discretization and (b) isotherm 
distribution of 5-subsection 3D rectangular slab. 

The plot of the L 2 norm progression is shown in Figure 19 for the first 10 iterations, however, 
the level of convergence of e ■ AT was achieved after only 4 iterations. Table 2 reveals the 
memory requirement proportions for each case and the computation time for the algebraic setup 
and solution after the 4 iterations for which L 2 = 0.012. 



Figure 19. L 2 norm progression for the 5-region 3D slab. 
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Table 2. Memory and time comparison for 3D slab problem. 



1 -Region 

5 -Region 

Iterations 

0 

4 

L 2 Norm 


0.012 

Memory 


BEM 

Time for Setup 

488s 

16s 

Time for Solution 

374s 

86s 


5.2.3 Steady Conduction Model of a 3-Region 3D Thrust Vector-Control Vane 

A 3D Thrust Vector-Control Vane is initially discretized in a single-region using 610 
equally spaced bilinear isoparametric discontinuous boundary elements with a total of 2440 
nodes around the boundary. Figure 12. a. shows the BEM discretization for this problem along 
with the isotherm distribution in Figure 12.b. The boundary conditions were distributed as 
insulated on the bottom surface, convective on the back surface with = 200 and h = 100, 
and convective with T \ = 4000°(7 and h(x ) = 1000(x/x mas ) 2 W /m?K for a maximum of 
h = 1000 W /m 2 K on the leading edge. A thermal conductivity of k = 14.9 W /mKwas used. 



(a) 


(b) 


Figure 20. (a) BEM discretization and (b) Isotherm distribution of single-region 3D Vane. 
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Next, the 3D Thrust Vector-Control Vane is divided into 3 sub-domains and discretized with a 
maximum of 300 bilinear isoparametric discontinuous boundary elements with 1200 nodes. 
Figure 13. a. shows the BEM discretization for this problem and the isotherm distribution is 
plotted in Figure 13.b. for the non-converged solution after 1 iteration. Figure 14. a shows the 
non-converged solution isotherms after 40 iterations and finally Figure 14.b. shows the 


converged solution isotherms after 80 iterations. 



Figure 21. (a) BEM discretization and (b) Isotherm distribution of 3-region 3D vane after 1 


iteration. 



(a) 


(b) 


Figure 22. Isotherm distribution of 3-region 3D vane after (a) 43 iterations and (b) 81 iteration. 
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The plot of the L 2 norm progression is shown in Figure 23 comparing the iterations for the cases 
of adiabatic initial guess at the subsection interfacial nodes and physically-based initial guess. 
When an adiabatic intial guess is made at the sub-structure intefaces, it took 81 iterations to reach 
the level of convergence of e • AT, while the physically-based intial guess provided an intial 
error norm of 0.045 (vs 0.32 with the adiabatic guess) which lead to less than 45 iterations to 
achieve convergence. This clearly demonstrates the effectiveness of the proposed physically- 
based intial guess in reducing the computational load of the iterative process in this case. Table 3 
reveals the memory requirement proportions for each case and the computation time for the 
algebraic setup and solution. Although final time to solution is comparable, the larger problems 
could not be tackled by the 1 region approach due to memory requirements and round-off error. 



Figure 23. L 2 norm progression for the 3-region 3D vane comparing the iteration for the cases of 
adiabatic initial guess at the subsection interfacial nodes and physically-based initial guess . 

Table 3. Memory and time comparison for 3D slab problem. For the 3-region: (a) adiabatic guess 
and LU factors computed at each iteration, (b) adiabatic guess and LU factors computed once 
and stored, (c) Eq. (32) intial guess and LU factors computed once and stored. 



1 -Region 

3-Region (a) 

3-Region(b) 

3-Region(c) 

Iterations 

0 

81 

81 

43 

Memory 

100% 

24% 

24% 

24% 

Time for Setup 

605s 

391s 



Time for Solution 

916s 

15,324s 

793s 

421s 
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5.2.4 Modeling a 3-D Cooled Blade 

Here, a blade with a 10 cm chord and it is 14cm in the spanwise direction. The blade is 
cooled by two plena, see Fig. 24. The blade is discretized using GridPro™ [17] into 6 
subsections with a surface grid of a total of nearly 6, 000 bilinear elements or nearly 24, 000 
degrees of freedom, see Fig. 25. Each block is kept at a discretization level near to 1000 bilinear 
boundary elements. Adiabatic conditions are imposed on the top and bottom surfaces of the 
blade. Convective boundary conditions are imposed on all other surfaces. The film coefficient on 
the outer surface of the blade is taken as h = 1000 W /rrrK with the reference temperature 
taken as 1000 AT, while the cooling plena are both imposed with film coefficients h = 500 
W/m 2 K with the reference temperature taken as linearly varying from 300 A to 400 A" in the 
increasing ^-direction of the cooling plenum closest to the leading edge, while linearly varying 
from 500 A' to 400K in the decreasing ^-direction of the cooling plenum closest to the trailing 
edge. 
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Figure 24. BEM grid for 3D cooled blade. 


NASA/CR- 


- 2003-212195 


44 


Figure 25. Substructuring of a 3-D plenum-cooled turbine blade. 


The physically-based initial guess at the sub-section interfaces using Eqn. (32) provided 
an excellent starting point for the iteration which converged in 16 steps to provide an L 2 iterative 
norm, defined in Eqn. (41), of L 2 = 0.00011698. It took 34,905s to set up the matrices, obtain 
and store their LU factors, and 813s to solve the problem iteratively. In this case, the LU factors 
were stored for each sub-section BEM model, and each iteration consisted of a forward and 
backward substitution. The resulting temperature plots illustrated in Fig. 25 and Fig. 26 reveal a 
very smooth distribution across all blocks. The resulting surface heat fluxes are presented in Fig. 
27 revealing a very smooth distribution from a minimum of — 180, 000 W /m 2 K to a maximum 
of 230,000 W/m 2 K. It should be noted that the subsectioning approach is ideally suited for 
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parallel implementation. The authors are pursuing this avenue prior to integration of the 
algorithm with the CHT solver. This concludes the example section. 


Y 



T: 650 700 750 800 850 900 950 1000 


Figure 26. Converged surface temperature distribution ( K ) 
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X 



Q: -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 


Figure 27. Converged surface heat flux distribution (j^W/m 2 K ) 


6 Conclusions 

A combined BEM/FVM approach using the TFFB conjugate method has been 
implemented in a 3-D context to model CHT in cooled turbine blades. As a boundary-only grid is 
used by the BEM, the computational time for the heat conduction analysis is insignificant 
compared to the time used for the NS analysis. The proposed method produces realistic results 
without using arbitrary assumptions for the thermal condition at the conductor surface. Results 
from a CHT numerical simulation of a 3-D film-cooled blade section are presented and are 
compared with those obtained from the standard approach of a two temperature model. A 
significant difference in the level and distribution of the metal temperatures is found between the 
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two models. These differences have severe implications in materials design and subsequent 
thermal stress analysis of the blade carried out using these metal temperatures. In practice, 
turbomachinery components such as modern cooled turbine blades which often contain several 
hundred film cooling holes and intricate internal serpentine cooling passages with complex 
convective enhancement configurations such as turbulating trip strips. This poses a real 
computational challenge to BEM modeling. The subsectioning iterative approach outlined in this 
report thus offers a promising technique to address this problem. 

Future Work 

We are in the process of building a PC-based cluster running under Windows. We are 
extending the sub-structure code under MPI and MPI-2 parallel protocols. Testing and fine 
tuning of the code will be undertaken in the course of the next year. We will implement a 
GMRES iterative solver to speed up computations at the sub-section level of the current code. 
This work is expected to continue into January. We will also begin to develop a 3D parallel 
multipole BEM fast solver to be run with a GMRES pre-conditioned non-symmetric solver. This 
phase of the project will be undertaken over the next two years. 
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Appendix: Numerical Evaluation of the Influence Coefficients 


The process is most readily illustrated by considering constant elements. In the constant 
element, which is a subparametric element, the field variables, T and q, are modeled as constant 
across each element while the geometry is represented locally as bilinear planes. Figure A.l. 
below shows a typical constant boundary element along with its transformed representation in the 
local rj — C coordinate system. 



Figure A.l. Bilinear subparametric boundary element. 


Notice that the geometric nodal locations of the element are ordered counterclockwise 
such that the normal vector always points outwards from the domain of the problem. The global 
coordinate system (x. y, z) is transformed into a local coordinate system (//. (') using the bilinear 
shape functions as given in Eqn. (15). However, the temperature and heat flux are modeled as 
constant with the node located at the geometric center of the boundary element, thus 

T j(v,0= T j and Qj(v,0 = Qj (A.l) 
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Clearly, the temperature and heat flux are discontinuous at the element interfaces. Thus, a 
constant element is termed discontinuous. Introducing the above discretization in the BIE and 
noting that NPE = 1, and collocating the discretized BIE at each of the boundary nodes £, there 
results 


N N 

j = 1 J =1 


where the influence coefficients Hij and G tJ are defined as 


Hij = ® q*(xj,£i)dS(x ) 


AS, 


Gij = f T*(xj,£i)dS(x ) 


AS, 


(A. 2) 


(A. 3) 


These are evaluated numerically using Gauss-Legendre quadratures with an adaptive scheme to 
be discussed shortly. Although very simple in implementation, use of the above constant element 
formulation does not lead to satisfactory results in many cases, and these are discussed purely for 
illustration. Introducing the definition of bilinear representation of the geometry into the constant 
element influence coefficient definition, the constant element influence coefficient integrals are 
explicitly, 



C)cMC 

T*(x(ri , (),&)Jj{v,()dv d( 


(A.4) 


The Jacobian of the transformation over the j-th element, J ? (r/, (') is 

Jj(vX) = 9 m 9<( ~ 


(A. 5) 


where g m , , and g r/ > are the components of the metric tensor defined as, 
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9r}i] — 

9(( = 
9vC = 


r\ \ 2 / r\ \ 2 / r\ \ 2 

ox \ i oy\ ( oz \ 
dr) ) + \dr) ) + {dr)) 

2 



/ dx dx\ ( dy dy\ fdzdz\ 

\dr) d() + \drjdC) + \dr) d( ) 


(A. 6) 


The metrics §^, §|, and §4 are readily found by differentiation of Xj(rj, (), 

Uj(r 7 , C), and Zj(r), () in Eqn. (15). Introducing the metrics into the Jacobian and simplifying 
leads to, 


Jj(v, 0 = (A.7) 

I / dy dz dz 8y\ 2 7~dz~dx dx dz\ 2 7 ~dx~dy dy dx\ 2 

\j \dr) d( dr) d( ) \drj d( dr) d( ) ^ \dr) d( dr) d( ) 


In the expression for Q*(x, & there arises the need to evaluate the outward-drawn normal n, as 
by definition 

^ t dT*(x,^) )V7 *, ^ ^ 

9 (x, &) = - k = - kVq (x, &) ■ n (A.8) 

a n 


The components of the unit vector on each elements can be easily computed using, 


n x = 
n y = 
n z = 


( dy dz 
\dr) d( 
f dz dx 
\ dq d( 
f dx dy 
\dr) d( 


dz dy\ 
dr) d( ) 
dx dz\ 
dr) d() 
dy dx\ 
dr) d() 


A (v, 0 

A Oh 0 


(A. 9) 


The numerical integration process necessary to obtain the influence coefficients is performed by 
double Gaussian quadratures (Gauss-Legendre specifically) simultaneously along the two local 
axis r) and leading to the following form, 
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(A. 10) 


NG NG 


Ha = Y Y w - w ^i x jiVumt C nm)i £,i)Jj(jlnmi Cnm , ) 
n—1 m— 1 
NG NG 

G ij = Y Y WnWmT * ( x j ( Vnm, 5 Cnm, ) i Vnm . > Cnm , ) 


n=l m.= 1 


where NG is the number of Gaussian points or order of integration employed, rj nm and Cnm, are 
the locations of the Gaussian quadrature points (zeroes of the appropriate Legendre polynomials), 
and W n and W m are the quadrature weights [58,59]. 

The number of Gaussian points employed can be adapted to every integral depending on 
the variability of the integrand. The influence coefficients G are inversely proportional to the 
Euclidean distance between the field or integration point x and the collocation point Ci, and the 
influence factors H are inversely proportional to the square of the Euclidean distance between 
the field or integration point and the collocation point Ci- Therefore, as the collocation point Ci is 
positioned closer to the integration element, the variability of the integrand increases requiring an 
increase in the number of Gaussian points, NG, for the integral approximation to provide a 
similar level of accuracy. Hence, a simple distance rule can be heuristically employed to change 
the number of Gaussian points depending on how far the collocation point is to the integration 
element. For non-singular elements, the following ratio of lengths is computed prior to 
integration 


r c 


T ij—max 
G 'ij—min 


(A. 11) 


where r, y is measured from the collocation point to 49 (7x7 equidistant) locations along the 
integration element j. The following is a table of heuristic quadrature and adaption rules 
determined by experience that are adopted in the 3-D BEM code: 
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Table A.l Gauss and adaption rules used in the BEM code. 


r c Range 

# of Integration Cells 

Gaussian Points 

1 < r c < 1.2 

1 

9 (3 x 3) 

1.2 <r c < 1.5 

1 

25 (5 x 5) 

1.5 <r c < 2.0 

1 

49(7x7) 

2.0 < r c < 3.0 

1 

100 (10 x 10) 

3.0 < r c < 5.0 

1 

225 (15 x 15) 

5.0 < r c < 7.5 

1 

400 (20 x 20) 

7.5 <r c < 10 

1 

625 (25 x 25) 

10 < r c < 25 

4(2x2) 

4 x 900 (30 x 30) 

25 <r c < 50 

9 (3 x 3) 

9 x 900 (30 x 30) 

50 < r c < 100 

16(4x4) 

16 x900 (30 x 30) 

100 <r c < 200 

25 (5 x 5) 

25 x 900 (30 x 30) 

200 <r c < 500 

36 (6 x 6) 

36 x 900 (30 x 30) 

r c > 500 

64 (8 x 8) 

64 x 900 (30 x 30) 


However, simply increasing the number of Gaussian integration points is in some cases not 
enough to obtain an accurate approximation to the integral. This is precisely the case when the 
collocation point £, is extremely close to the integration element. In such cases a 
subsegmentation of the element is required and can be performed in a similar fashion as to the 
increase of Gaussian points, that is, the closer the collocation point is to the element, the more 
subdivisions are made to the element with a fixed number of Gaussian points for each 
subelement. 

The particular case in which the collocation point is located over the integration 
element must be treated with caution. For this case the integrand becomes singular lacking an 
accurate integral approximation through a regular Gaussian rule. Even though the integral for the 
influence factor Ha is strongly singular because its integrand is inversely proportional to the 
distance squared (1/r 2 ), this integral need not be computed directly unless regularized, and 
instead it can be evaluated using the equipotential relation by summing the off-diagonal terms 
[44-46], that is 


N AS A/CR— 2003 -212195 


59 












(A. 12) 



The integral for the influence coefficient Gu is however weakly singular because its integrand is 
inversely proportional to the distance (1/r), and can be accurately computed through a 
transformation of the local coordinate system which effectively clusters the Gauss points. Figure 
A. 2 shows the primary subsegmentation of the singular element. Twelve quadrilateral 
subdivisions have been made to the singular subparametric boundary element in Fig. A. 2. The 
shaded area "0 " corresponds to the singular portion of the element that is to be transformed into 


a local polar coordinate system. 



Figure A. 2. Subsegmentation of singular - bilinear subparametric boundary element. 
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Figure A. 3. Polar coordinate transformation of the area "0" of the singular boundary element. 


The shaded area of the singular element is fit into a new local system ( 77 , () and subdivided into 
four triangular subelements. Each of the subelements is transformed into a new local polar 
coordinate system (p, 9) where, 

77 = p cos 9 (A. 13) 

( = p sin 9 


therefore, the factor which corresponds to the integral over the shaded area "0" is comprised 
of four integrals (G' ? 0 /, Gf, G° 3 , and G° 4 ) as, 


/^ y 01 


r> 02 


^ y 03 


04 

^ ii 



jJO 
.3 2 

4 / sin 


^ 57T 1 

4 C cosO 


do 


T*{xi(rj, C C)pdpd9 
T*(xi(v, 0 ,0)Mv, Op dp d9 
T*(xi{r 7 , C), CWp d9 



T*{xi(r), C ),£i)Mv, 0 pdpd9 


(A. 14) 
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Notice that the transformation of the differential area introduced the variable p to the integrand 
which relaxes the singularity of the fundamental solution T*(xj(rj, £),£,). In addition, an extra 
transformation is necessary to fit the limits of integration into the ( — 1,1) range such that. 


Pi 

P‘2 

P3 

P4 


S + 1 

2cos 9\ 
s + 1 

2 sin 02 
■s + 1 

2cos #3 
s + 1 

2 sin 04 


, 02 — -^(t + 2 ) 
, *3 = ^ + 4) 
, 0 a = -^(t + 6) 


(A. 15) 


where the subindeces of the coordinates p and 9 correspond to the transformation for the 
particular integral term GG , GG, and G®f. Therefore, the influence coefficient subintegrals 
are transformed into, 



T*(xi(rp C ),£i)Mv, 0 
T*(xi(rp C ),£i)Mv, 0 
T*(xi(rj, C) C) 
T*(xi(r ) , C ),£i)Ji(v, 0 


?rpi 

8cos 9 \ 
np2 

8sin 02 
^pz 

8 COS 03 

npA 

8sin 04 


ds dt 
ds dt 
ds dt 
ds dt 


(A. 16) 


which can be solved numerically by the use of Gaussian integration and through the 
corresponding transformations detailed in the previous relations. Each triangle is integrated using 
a 100 (10 x 10) Gaussian point rule distributed in the transformed r — 9 system. Each of the 
remaining 12 cells that do not contain singularities are integrated using a 900 (30 x 30) Guassian 
point rule. 


When dealing with isoparametric bilinear elements, the procedure of evaluation of the 
element influence coefficients Hf- and G ■ is the same, except that the shape functions M k (>],(), 
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k = 1...4, appear multiplying each of the integrands of the influence coefficients. That is for 
instance, 

G ij = J' J_T’(x 1 ( V X)-X,)M t ( V ,i)J 1 (vX)dvdi; (A. 17) 

NG NG 

G{j = ^ ^ W n \V m T (xj(ri nm i C, n m) i ' (j] nm: C, nrn ')Jj(v n ! mi Cnm,) 

n—1 m— 1 

where the subscript i = L..4jY refers to the collocation point which is positioned at all four 
offset nodal locations of each element. The subscript j = 1...N refers to the integration element 
and the subscript k = 1...4 corresponds to each of the shape functions M k (r)£) in the integrand. 

The subsegmentation described for the G %% influence coefficient in the subparametric 
element is extended to the case of discontinuous isoparametric elements by performing a non- 
symmetric subsegmentation of the element i depending on the location of the collocation point 
on the boundary element, see Fig. (A. 4). 
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(•) Collocation T&q node 
O T&q node 
• Geometry node 


Figure A. 4. Subsegmentation of singular bilinear isoparametric boundary element. 
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